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We show that in large-field inflationary scenarios, superheavy (many orders 
of magnitude larger than the weak scale) dark matter will be produced in cos- 
mologically interesting quantities if superheavy stable particles exist in the 
mass spectrum. We show that these particles may be produced naturally dur- 
ing the transition from the inflationary phase to either a matter-dominated 
or radiation-dominated phase as a result of the expansion of the background 
spacetime acting on vacuum quantum fluctuations of the dark matter field. 
We find that as long as there are stable particles whose mass is of the order 
of the inflaton mass (presumably around 10 13 GeV), they will be produced in 
sufficient abundance to give fi = 1 quite independently of any details of the 
non-gravitational interactions of the dark-matter field. 
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I. INTRODUCTION 



It is now commonly accepted that most of the mass in galactic halos as well as in the 
Universe as a whole is composed of dark matter (DM). There are many indications that 
the DM consists of some new, and yet undiscovered, weakly interacting massive particles 
(WIMPs). 

Despite the fact that the nature of the DM is still unknown, it is usually thought that 
DM particles cannot be too heavy. If the WIMP is a thermal relic, then it was once in 
local thermodynamic equilibrium (LTE) in the early universe, and its present abundance 
is determined by its self-annihilation cross section. From unitarity arguments one 
expects the mass of a thermal relic to be less than about 500 TeV. The present abundance 
of non-thermal relics is not determined by their self-annihilation cross section since they 
needn't have been ever in LTE in the early universe. An example of a non-thermal relic 
is the axion, and the present axion abundance is determined by the dynamics of the 
phase transition associated with symmetry breaking. Non-thermal relics are typically 
very light, e.g., the axion mass is expected to be in the range 10~ 5 to 10 _2 eV 0. 

Because the assumption of relatively low-mass DM seems quite natural, it is rarely 
questioned|] The goal of this paper is to show that the Universe might be made of 
superheavy WIMPs (we will refer to them as X particles), with mass larger than the weak 
scale by several (perhaps many) orders of magnitude. Two conditions are necessary for 
this to happen: a) the X particles must be cosmologically stable and b) their interaction 
rate must be sufficiently weak such that thermal equilibrium with the primordial plasma 

was never obtained. This second condition is easy to satisfy as long as the particle is 

Of course, superheavy dark matter particles have been considered before to a certain extent. In 
particular, there is an extensive literature regarding observational constraints on unusually heavy dark 
matter candidates (for example, see Refs. [[|, [^), and references therein). However, they do not 
restrict our scenario nor do they consider our production mechanism. 
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extremely massive (of the order of the Hubble parameter at the end of inflation) . 

We point out that superheavy dark matter may be created during the evolution of 
the Universe in a number of ways. If it is produced during the process of reheating after 
inflation, then the upper bound on its mass Mx can be as large as the reheating temper- 
ature Trh- The latter should be less than about 10 9 GeV in order to avoid overproducing 
dangerous relics such as quasistable gravitinos in supergravity inspired scenarios. The 
mass upper bound can be pushed higher than the reheating temperature if one allows 
the DM to be produced directly through the decay of the inflaton field. In that case, 
the mass upper bound is the inflaton field mass, which is presumably less than about 
10 13 GeV. On the other hand, if reheating after inflation is preceded by a preheating 
stage || it is certainly possible to produce by resonance effects copious amounts of dark 
matter particles with masses much larger than the inflaton mass . 

In this paper, we consider yet another mechanism of generating heavy DM. We study 
the possibility that DM is produced in the transition between an inflationary and a 
matter-dominated (or radiation-dominated) universe due to the "nonadiabatic" expan- 
sion of the background spacetime during the transition acting on the vacuum quantum 
fluctuations. 

The distinguishing feature of this mechanism is the capability of generating particles 
with mass of the order of the inflaton mass (usually much larger than the reheating 
temperature) even when the particles only interact extremely weakly (or not at all) with 
other particles and do not couple to the inflaton(s). We find that they may still be 
produced in sufficient abundance to achieve critical density today due to the classical 
gravitational effect on the vacuum state at the end of inflation. More specifically, we 
will show that in the range 0.04 < Mx/H < 2, where H ~ ~ 10 13 GeV is the 
Hubble constant at the end of inflation (m^ being the mass of the inflaton), the DM 
produced gravitationally can have a density today of the order of the critical density. 
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This result is quite robust with respect to the "fine" details of the transition between 
the inflationary phase and the matter-dominated phase, and independent of the coupling 
of the DM to any other particle. This result is reasonably robust also with respect to 
the ambiguity associated with the choice of the vacua as we have tried to minimize the 
number of particles produced by choosing an infinite adiabatic order in-out vacua. The 
only "non-trivial" requirements, other than that large field inflation occurs, are that the 
WIMPs posses a mass close to the inflaton mass and that they are stable. 

Mechanically, the DM particle creation scenario is similar to the inflationary genera- 
tion of gravitational perturbations that seed the formation of large scale structures (see 
for example the review given in Ref. ||). In the usual scenarios of this form, however, 
the quantum generation of energy density fluctuations from inflation is associated with 
the inflaton field which dominated the mass density of the universe, and not a generic, 
sub-dominant scalar field. 

Because it is usually assumed that DM forms from the decays or interactions of the 
reheating products, it usually has a stage of LTE in its early history. In our scenario 
the large mass of the dark-matter particle will prevent it from thermalizing, and its 
abundance will depend only on its mass and the behavior of the spacetime, not on its 
weak coupling to other nongravitational fields. 

Others have considered gravitational particle production at the end of inflation. For 
example, Ford [(J and Yajnik [|10J both consider particle production as a result of the 



nonadiabaticity of the transition from an inflationary phase to a matter or radiation dom- 
inated phase (although with a different cosmological implication in mind). Ford treats 
only massless, non-conformally coupled fields using a well known perturbation technique 
(see references within IJ). Yajnik considers minimally coupled scalar field theory in the 
limit of small masses, with an abrupt transition from an inflationary phase to a radiation 
dominated phase. In our work, we consider extremely massive, conformally coupled fields 



and calculate the particle production exactly by numerically solving the mode equation. 
We treat the conformally coupled case because conformal coupling generally minimizes 
the number of particles produced, particularly in the small mass ranges. Unlike Yajnik, 
we also consider the case where the metric is an analytic function of the conformal time 
and show that this leads to qualitatively different behavior of the density of particles 
produced for large masses. The analyticity implies a conservative estimate since fewer 
particles are produced in that case than in the abrupt transition case. 

Some of the ideas present in our scenario are also contained in the work of Linde and 
Kofman [11], [12], and [|H|. However, the purpose of their work was to point out that 
isocurvature cosmological (large scale) perturbations can be produced during inflation. 
They did not consider the importance of the nonadiabaticity of the transition at the 
end of inflation which is responsible for the production of our superheavy dark matter. 
Instead, they mainly relied upon estimates of the particle production during the de Sitter 
phase or the classical (long wavelength) component of the particle field energy density 
left over after inflation. 

This paper is organized as follows. In the next section, we elaborate on the dark 
matter scenario and the calculational method. In Section III, we discuss the numerical 
results. We then summarize our work in Section IV. In the appendix, we derive the 
asymptotic mass dependence of the dark matter density presented in Section II. 



II. SCENARIO AND CALCULATIONAL METHOD 



In this section we discuss the dark matter abundance calculation in our scenario. 
First, we give an expression for the dark matter density today in terms of the number 
density when it was produced. We then consider the mass range of the dark matter 
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necessary if it is never to thermalize. Finally, we discuss the mechanics of the gravita- 
tional production of particles. In particular, we discuss the number density definition 
and present the asymptotic dependence of the number density on the particle mass. 

Suppose the dark matter never attains LTE and is nonrelativistic at the time of 
production. The usual quantity Qxh 2 associated with the dark matter density today can 
be related to the dark matter density when it was produced. To develop the relation, we 
begin by writing 

Pn(to) ~ Pr^rh) V T ) ' U 
where denotes the energy density stored in radiation, p x denotes the energy density 
residing in the dark matter, Trh is the reheating temperature, To is the temperature 
today, to denotes the time today, and tun denotes the approximate time of reheating 
completion.^ To obtain pxitiai) I Pr^bm), we must determine when X particles are 
produced with respect to the completion of reheating and the effective equation of state 
operative between X production and the completion of reheating. 

At the end of inflation the universe may have a brief period of matter domination 
resulting either from the coherent oscillations phase of the inflaton condensate or from 
the preheating phase |J. If the X particles are produced at time t = t e when the de 
Sitter phase ends and the coherent oscillation period just begins, then both the X particle 
energy density and the inflaton energy density will redshift at approximately the same 
rate until reheating is completed and radiation domination begins. Hence, the ratio of 
energy densities preserved in this way until the time of radiation domination is 

p R (t RH ) ~ 3 M Pl 2 H 2 (t e y { ' 

where Mpi ~ 10 19 GeV is the Planck mass and most of the energy density in the universe 

just before time tRH is presumed to turn into radiation. Thus, using Eq. ([[]), we may get 
6 More specifically, this is approximately the time at which the Universe becomes radiation dominated. 
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an expression for the quantity Qx = pxifo)/ Pc{to), where pc(to) = 3HqMp 2 /8n and 
H = 100 h km sec" 1 Mpc" 1 : 



Here VL^h 2 « 4.31 x 10~ 5 is the fraction of critical energy density that is in radiation 
today and nx is the density of X particles at the time when they were produced. 

Note that because the reheating temperature must be much greater than the temper- 
ature today (Trh/Tq > 4.2 x 10 14 ), in order to satisfy the cosmological bound Q x h 2 < 1, 
the fraction of total energy density in the dark matter at the time when they were pro- 
duced must be extremely small. To illustrate this, take H 2 {t e ) ~ m 2 , and p(t e ) ~ m^Mp 2 . 
Then Q x h 2 ~ 10 17 (T RH / 10 9 GeV)(p x (t e ) / p(t e )) . It is indeed a very small fraction of the 
total energy density we wish to extract in the form of massive X particles. 

This means that if the dark matter particle is extremely massive, the challenge lies 
in creating very few of them naturally. We will see that the gravitational production 
naturally gives the needed suppression. Note that if reheating occurs abruptly at the end 
of inflation, then the matter domination phase may be negligibly short and the radiation 
domination phase may follow immediately after the end of inflation. However, this does 
not change Eq. @. 

For the superheavy X particles to be good candidates for DM, they have to be stable 
or at least have a lifetime greater than the age of the universe. This may occur in 
super symmetric theories where the breaking of supersymmetry is communicated to the 
ordinary sparticles via the usual gauge forces flip . In gauge- mediated supersymmetric 
models there are two sectors with possible stable particles which might act as superheavy 
dark matter candidates: 

1) The secluded sector, which is strongly interacting: Supersymmetry is broken dy- 
namically and some F-term gets a nonvanishing expectation value, where the scale of 




(3) 
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super symmetry breaking, as usual, is denoted by yF. 

2) The messenger sector: This sector contains the fields charged under the SU(3)c <S> 
SU{2)l ® U{1)y gauge interactions, and communicate supersymmetry breaking to the 
sparticles in the observable sector. The mass of the messenger fields is usually denoted 
by M. 

After the messengers have been integrated out, sfermions receive a mass squared 
m 2 ~ a 2 A 2 , where a is the appropriate gauge coupling and A ~ F/M. Notice, in 
particular, that the spectrum of the superparticles depends on the ratio A = F/M which 
is fixed to be relatively small and in the range 10 to 10 3 TeV. However, this does not 



necessarily mean that yF and M are of the same order of magnitude as A [15 since it 
is only their ratio which is fixed around 10 3 TeV: the hierarchy yF, M 3> A is certainly 
allowed ]16 . 



The secluded sector often has accidental symmetries analogous to the baryon number. 
This means that the lightest particle in the secluded sector might be stable and a good 
candidate for dark matter with a mass of the order of yF, much larger than the weak 
scale. The lightest messenger field might also be a good candidate for superheavy DM. 
Indeed, if the supersymmetry breaking sector contains only singlets under the SU(3)c ® 
SU{2)l (gi U(l)y gauge interactions and if there are no direct couplings between the 
ordinary and messenger sectors, then the theory is characterized by a conserved global 
quantum number carried only by the messenger fields. The typical mass M of the DM 
component in the messenger sector may be much larger than the weak scale. 

Another framework in which we might expect the presence of superheavy stable 
particles is a Kaluza-Klein theory (a unified theory which requires space-time dimensions 
higher than four). A popular example is provided by M-theory |I7J where the number 
of dimensions is D = 11. These theories are characterized by the presence of a tower 
of Kaluza-Klein modes which are left after the compactification of the extra D — 4 



dimensions. For instance, if D = 5, the existence of a compact fifth dimension implies 
an infinite tower of four- dimensional particles corresponding to quantized excitations of 



the extra dimension. These massive particles have been called "pyrgons" If any 

of the pyrgon states are stable or have a lifetime greater than the age of the universe, 
they might act as DM with a mass of the order of the inverse of the physical size of the 
compact dimensions -R^ 1 , which is likely to be larger than the weak scale by many orders 
of magnitude. 

For the gravitational production scenario to be distinguishable from other scenarios, 
X must never thermalize. The condition for the dark matter particles to be out of 
equilibrium and their comoving number density to be constant is 

n x (a A \v\)<H, (4) 

where H is the Hubble parameter, (a a \v\) is the thermal averaged self-annihilation cross 
section times the M0ller speed for the dark matter particles X. Since the cross section 
a a is expected to be at most about M^ 2 (usually smaller; sometimes much smallerf]) and 
nx is bounded by the condition that Qxh 2 < 1, we obtain from Eq. (^) 

nxMvD _ 7 x IP' 19 (g/MpQ 

H ~ (7W/10 9 GeV) (M X /M P1 )3 1 } 

as the quantity which must be less than one at t — t e to avoid thermalization. For a low 
reheating temperature of 10 2 GeV and a typical value of H = 10~ 6 Mpi for inflationary 
scenarios, we find a conservative condition Mx/H > 1 for the particles never to reach 
chemical thermal equilibrium. Note that this is a rather conservative estimate since the 
reheating temperature is likely to be larger and the cross section is likely to be smaller. 

We also remark that because the reheating temperature is likely to be much smaller than 

7 For example, if there is a heavy gauge particle mediating the process, then the effective coupling 
will be further suppressed and the relevant mass scale for the cross section will be the mediating particle 
mass instead of the X mass. 
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the X mass, the thermal production of the X particles is negligible.^ 

Now let us describe the basic physics underlying our mechanism of gravitational 
production of DM. 

In this paper we take space-time both in and out of inflationary era to be spatially 
flat, homogeneous, and isotropic, with the line element of the form 

ds 2 = a\r 1 )(drf -rfx 2 ). (6) 

For simplicity (and without much loss of generality), we restrict ourselves to a massive 
scalar field coupled to classical gravity and nothing else. The other couplings are assumed 
to play an insignificant role in the gravitational production. 

There are various inequivalent ways of calculating the particle production due to 
interaction of a classical gravitational field with the vacuum (see for example ||19|| , pO 



and |21[). In our work, we use the method of finding the Bogoliubov coefficient for the 
transformation between positive frequency modes defined at two different times. We will 
show below that the large mass dependence of the DM number density is determined 
by either the differentiability (or the smoothness) of the scale factor or the choice of 
the vacuum. On the other hand, for Mx/H < 1 where H is the value at the end of 
inflation, the results are quite insensitive to the differentiability or the fine details of the 
scale factor's time dependence. For 0.04 < Mx/H < 2, we find that all the dark matter 
needed for closure of the universe can be made gravitationally, quite independently of 
the details of the transition between the inflationary phase and the matter dominated 
phase. 

To see the effects of vacuum choice and the scale factor differentiability on the large 
X mass behavior of the X density produced, we start with the canonical quantization 



8 Since for times larger than t e , the interaction rate continues to be smaller than H, the particles will 
not thcrmalizc later cither. 
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of the X field in an action of the form (in the coordinate ds 2 = dt 2 — a 2 (t)dx 2 ) 
« 3 (<r2 (VX) 2 



s 



J dt J d ^^.lx 2 - - M 2 X X 2 - £RX 2 ) (7) 



where R is the Ricci scalar. After transforming to conformal time coordinate, we use the 
mode expansion 

d 3 k 



(27r) 3 / 2 a(r/) 



where because the creation and annihilation operators obey the commutator [a^, aj.J = 
<5®(ki — k 2 ), the h^s obey a normalization condition hkh' k * — h' k h* k = i to satisfy the 
canonical field commutators (henceforth, all primes on functions of r\ refer to derivatives 
with respect to if). The resulting mode equation is 

K(v)+™l(v)h k (v) = 0, (9) 

where 

w 2 k = k 2 + M\a 2 + (6f - l)a"/a . (10) 

The parameter £ is 1/6 for conformal coupling and for minimal coupling. From now 
on, we will set £ = 1/6 for simplicity but without much loss of generality. By a change 
in variable r\ — > k/a, one can rewrite the differential equation such that it depends only 
on H(r)), H'(r))/k, k/a(rj), and Hence, we introduce the parameter Hi and a« 

corresponding to the Hubble parameter and the scale factor evaluated at an arbitrary 
conformal time rji, which we take to be the approximate time at which X are produced 
(i.e., rji = T)(t e )). We then rewrite Eq. @ as 

h'{{fj) + (P + & 2 5 2 )/^) =0 [b = Mx/H,} (11) 

where rj = rjciiHi, a = aja^ and k = k/{fliHi). For simplicity of notation, we shall drop 

all the tildes from now on. This differential equation can be solved once the boundary 
9 This differential equation is h%(y) + (l/H)H'(y)h' k (y) + (1 + M x /y 2 )/H 2 (y)h k = 0, where y = k/a. 
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conditions are supplied. Since the annihilation operator is just a coefficient of an ex- 
pansion in a particular basis, fixing the boundary conditions is equivalent to fixing the 
vacuum. 

To obtain the number density of the particles produced, we will perform a Bogoliubov 
transformation from the vacuum mode solution with the boundary condition at rj = 770 
(the initial time at which the vacuum of the universe is determined) into the one with the 
boundary condition at i) = T71 (any later time at which the particles are no longer being 
created). In the examples given in the next section, 770 will be taken to be — 00 while r\\ 
will be taken to be at +00 in order to define vacua of infinite adiabatic order (explained 
below) which results in a smaller particle production than for any finite adiabatic order 
vacua.0 The exact values of 770 and rji are not important for those examples as long 
as they are in a region in which a' /a 2 <C 1 or ba/k <C 1. Defining the Bogoliubov 
transformation as h^irj) = akfi^(rf) + /3kh* k Vo (if) (the superscripts denote where the 
boundary condition is set), we have the following energy density in the particles produced: 

/ 1 \ 3 poo au „ 

px(vi) = M x n x ( Vl ) = M X H? l^—j jf ^k 2 \(3~f, (12) 

where^] one should note that the number operator is defined at 771 while the quantum 
state (approximated to be the vacuum state) defined at rjo does not change in time in 
the Heisenberg representation. 

As usual, there is an ambiguity in the definition of the vacuum, which is equivalent 
to an ambiguity in the boundary conditions of Eq. ([|). One method of systematically 
classifying the various inequivalent vacuum states is through the adiabatic vacuum |22| 
definition. The adiabatic vacuum definition allows one to construct and classify a set 

of mode equation solutions which reduce to the usual plane waves when 0/(77) = for 

10 In the numerical calculation, one can only approximate these infinities with large numbers, but the 
limit is not singular. 

11 Here we restored the tildes for clarity. 
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all 77. The classification is based on a type of WKB asymptotic expansion in powers of 
conformal time derivatives of Wk- In particular, the classification allows one to quantify 
how two solutions with different boundary conditions (hence two vacua) will differ in 
terms of derivatives of w^,. Each derivative with respect to the conformal time is assigned 
a bookkeeping small parameter, and this small parameter's power in an expansion is 
referred to as the adiabatic order. We define AtYt adiabatic (order) vacuum at time rf 
by using the boundary condition 

= K(v*) = h'i A \ V *), (13) 

where (v) is a systematically chosen approximate solution to the mode equation that 
satisfies the mode equation up to Ath adiabatic order in the asymptotic limit that the 
adiabatic parameter goes to zero. Roughly speaking, the larger the adiabatic order of 
the vacuum, the closer it is to the Minkowski vacuum in the sense that it is less (in the 
adiabatic limit) dependent on the time at which it is defined. We refer the reader to the 
appendix (or Ref. |20]) for a more precise definition. 

As shown in the appendix, the asymptotic behavior of the number density as b —>■ oo 
can be obtained by the following rule: If the vacuum at r/o corresponds to an nth adiabatic 
vacuum, and the vacuum at 7/1 corresponds to a pth adiabatic vacuum, then as b — > 00 
the number density will behave likef^] 

n x ~ 6"( 2r+1 ) (14) 

where r = MinQo, n) provided that (d u a/dr] u )/a u+1 < 00 for all 77 G [r/ , ?7i] and all natural 
numbers z/.0 It is important to note that for a fixed time, the asymptotic expansion 

generated by Eq. ( |A2|) (in the Appendix) generally only converges up to a finite order. 

12 This behavior is also noted on pg. 69 of [^0) although there it is arrived at differently than in our 
Appendix. 

13 Note that by definition given in the appendix, p and n can be only even natural numbers or 0. 
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Hence, except under special circumstances an adiabatic vacuum of only a finite order can 
be generated. This means that in general, the number density will fall off with a finite 
power of 1/6 for large b. Only when an adiabatic vacuum of infinite adiabatic order can 
be generated, which usually means that the domain of a(t]) can be extended to ±00 with 
the property given above, does the number of particles produced fall off faster than any 
finite power of 1/6 (e.g., exponential suppression). In practice, we find that a spacetime 
which admits an infinite adiabatic order vacuum has the "advantage" of all the vacua 
defined in a sufficiently adiabatic region being numerically equivalent regardless of the 
vacua's adiabatic order and the exact time at which the vacua are defined. 

If within the domain there is one discontinuity of the first kincQ in (d g a/dr] q )/a g+1 
for q = s where —2 < s — 2 < r = Minfjo, n) and there are no discontinuities for q < s, 

n x ~b- {2s - z) . (15) 

This is true provided that (d v aj dr] u ) / a u+1 < 00 for all 77 in each of the continuous domain 
and all natural numbers v. Note that fractional power dependence on 1/6 will be possible 
if the discontinuity is not of the first kind (e.g., a(rj) = rf at rj = 0). 



III. NUMERICAL RESULTS 



We shall employ the method elaborated in the previous section to calculate the grav- 
itational production of particles in a couple of toy models of inflation which reflect some 
extreme ranges of differentiability. We will consider the case when the scale factor has a 
discontinuity of the first kind in one of its derivatives and when it is a C°° function. We 

will see that enough dark matter may be produced through this mechanism as to give 

14 The discontinuity of the first kind refers to the situation where the left and the right hand limits 
exist but are unequal. 
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critical density of dark matter today. 

Our example of a discontinuous model has the scale factor of the de Sitter space for 
i] < r\i and the matter or radiation dominated universe for 77 > 77^ 

M = l a i /(2- v /r h ) i£ V <Vi, m 
K ' \ ai(ri/rn) p if 77 > 77*, v ' 

where p = 2 for the matter dominated case and p = 1 for the radiation dominated case. If 
we define our vacuum states at t]q = —00 and rji = 00, then in this space time, adiabatic 
vacua of any order will be equivalent to infinite-order adiabatic vacua. For the transition 
into matter domination, the smallest derivative order in which there is a discontinuity 
(of the first kind) comes from a'(rj) at 77 — r\i while for the transition into the radiation 
domination, the analogous contribution comes from a"(rj) which has a discontinuity of 
the first kind at 77 — Hence, our analysis would predict that for large b, nx will 

increase like b in the case of matter domination whereas it will fall off like 1/6 in the case 
of radiation domination. 

Our second toy model looks at the other extreme limit of having a C°° function for 
a(rj) which behaves in the asymptotic limits of rj/r/i — > ±00 identically as the discontin- 
uous model: 

\ 1 - exp[-(r//^) 2 ] ( \ - tanh(r//^/2 + 1) \ ( V / Vl ) 2 ? 
a[V> ~ Q U (V/Vi) 2 { 2 y ) + ( i + 3exp( _^ ) )2 

. 1/2 

+ tanh(?7/^ - A) - tanh^/r?.; - 2A) > (17) 

where A = 1.07 is needed for proper normalization. The functional form was chosen 
rather arbitrarily except for the requirements of monotonicity, (d u a / dr] u ) / a u+1 < 00 for 
all 77 and natural numbers u, and appropriate power law asymptotic behavior. As before, 

this spacetime admits a vacuum of infinite adiabatic order at 77 = ±00. In Fig. [I] we show 

15 These discontinuities are unphysical and correspond at best to crude approximations. We consider 
them to test the sensitivity of our results on the choice of the model. 
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Figure 1: The scale factor (normalized to its value at r\ — T]i) is plotted as a function 
of the scaled conformal time rj/rji. The curves labeled "smooth" and "discontinuous" 
correspond to the toy models for the evolution of the scale factor without and with 
discontinuous derivatives, respectively. As rj/r/i — > — oo, the scale factor behaves as 
l/{rj/rji) 2 corresponding to a de Sitter space, and as rj/rji — ► +00, the scale factor behaves 
as (f]/i]i) 2 corresponding to a matter dominated Universe. The solid curve shows an 
analogous numerical solution of the scale factor for a (l/2)m 2 b (p 2 inflaton potential. 

how these models compare with the numerical result obtained by solving the (l/2)m^<j) 2 
inflationary model's equations of motion. Note that they differ mainly in the transition 
region (or the "nonadiabatic" region near f]/f]i = 1) where most of the particle production 
"occurs." 

In Fig. ^ we show the number density obtained numerically in these toy models. The 
peak at Mx/Hi ~ 1 for the C°° model is similar to the case presented in Ref. PBfl . As we 
expect, the large b behavior is exactly determined by the choice of the vacuum and the 
differentiability of the potential. Specifically, just as our asymptotic analysis showed, for 
large b, rix varies as b and 1/6 for the matter and radiation domination case respectively 
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Figure 2: The dark matter abundance today is shown as a function of the particle mass 
for various models. The mass is given in terms of Hi « 10~ 6 Mpi (the Hubble parameter 
at i] = rji, the beginning of the coherent oscillation period). In the "discontinuously 
into radiation" (if) has a discontinuity at r\ — rji, while in the "discontinuously 

into matter" case, a'(rj) has a discontinuity at rj — rji. The curves labeled "smoothly 
into" is for a(r]) that satisfies (d u a/drj u )/a u+1 < oo for all rj and natural numbers v. The 
curve labeled T = Hi/iOnz) shows a thermal density with this temperature. The un- 
shaded region satisfies the conservative nonthermalization condition obtained by setting 
{(J A \v\) = 1/M| in Eq. (|). 

in the discontinuous model, while nx is exponentially suppressed in the C°° model. 
Note that independently of the differentiability of the scale factor, if Mx ~ H% ~ 
for Trh ~ 10 9 GeV, X will have critical energy density today. On the other hand, 
unless there is some discontinuity in the scale factor for some nth derivative where n is 
not too large, this gravitational production mechanism will not generate enough dark 
matter in the Universe to give critical density for much larger masses (Mx 3> m$) even 
if such stable heavy particles exist. Furthermore, even in the mass range in which the 
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density of particles produced peaks, if the reheating temperature is below about 10 7 
GeV, this mechanism will most likely not generate significant amount of dark matter. 
In the case that this mechanism cannot produces these heavy particles, however, if these 
heavy particles couple to the inflaton and if preheating occurs, enough of them may be 
produced through the broad resonance mechanism to have critical density of superheavy 
dark matter today. 



IV. SUMMARY 



To conclude, we have investigated the scenario of creating nonthermalizing dark mat- 
ter gravitationally at the end of inflation (or the beginning of the coherent oscillation 
phase). There is a significant mass range ( O.lm^ to m^, where ~ 10 13 GeV) for 
which the X particles will have critical density today regardless of the fine details of the 
inflation-matter/radiation transition. Because this production mechanism is inherent in 
the dynamics between the classical gravitational field and a quantum field, it needs no 
fine tuning of field couplings or any coupling to the inflaton field. However, only if the 
particles are stable (or sufficiently long lived) will these particles give contribution of the 
order of critical density. For even larger dark matter masses, the broad resonance mech- 
anism of preheating (if it occurs) will produce these particles in sufficient abundance as 
to achieve Q — 1. 
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APPENDIX 

In this appendix, we derive Eq. (|i4|) and Eq. (|T5|), the asymptotic dependence of the 
dark matter density on the mass parameter b = Mx/Hi as b —>■ oo. This asymptotic 
behavior is, in general, dependent upon the choice of the vacuum state and the differen- 
tiability of the scale factor in an FRW type spacetime. We employ the adiabatic vacua 
ansatz p2| to classify the various possible (restricted) choices of vacua. Strictly speaking, 
our conditions for the various asymptotic behaviors are only sufficient conditions, but 
they have wide applicability as we demonstrate in this paper. 

Let us first review the concept of an adiabatic vacuum (see for example pg. 66 of 
Ref. [[2(|). We first define the concept of an adiabatic order as the power of 1/T that 
results for any term in a 1/T expansion after one makes the transformation r\ — > r\ and 
d/dt] — > T~ 1 d/di]. Note that if T _1 — > 0, then this is equivalent to an expansion in 
"smallness" of conformal time derivatives. The basic idea is that if the derivatives of the 
mode frequency are indeed small, then the degree to which the field theory breaks 
time translational symmetry can be characterized by the adiabatic order. This breaking 
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of the time translational symmetry^ is what is responsible for particle creation in our 
isotropic expanding Universe. 

To define the adiabatic vacuum, we first make a change in variables from h k to W k 
by writing 

1 



and obtain a new differential equation 



Wi 



wt - 1/2 \wz/w k -{3/2){wyw k 



(Al) 



(A2) 



where we have used Eq. (|TT1) and defined w k to be the coefficient of h~ k in Eq. (|lll)J I7 | 
Hence, let us define a map 



mi 



Ml 



\ 



2 1 








W k ~ ^ 




-1 


UtvJ 



(A3) 



which is a map that raises the adiabatic order by two and also define 



W, 



(n+2) 



A\W t 



(nh 



(A4) 



where the superscript denotes the adiabatic order and = w k . We can now write an 



approximate mode equation solution^ good to Ath adiabatic order as 



h 



(A) 



2W, 



(A) 



cxp 



-i fw^WW) ■ 



(A5) 



Finally, we define the adiabatic vacuum of Ath order at some value of r\ which we call 
rf by using the boundary conditionP] 



hM) = h[ A \r)% K( V *) = h' k (A \ V *), 



(A6) 



16 Conformal time translation generates conformal transformation in an FRW universe, and the mass 
term breaks the conformal symmetry. 

17 We have dropped all the tildes for simplicity in notation. Note also that a constant factor normal- 
ization choice of Eq. (Al) is unimportant for the Bogoliubov transformation. 

18 After finishing our paper, we learned that a complete and more precise analysis of the asymptotic 
behavior of the adiabatic modes can be found in Ref. p^ j. 

19 In the spirit of the adiabatic expansion, the equality needs to only be enforced to ^4th adiabatic 
order. However, we will for simplicity of argument assume throughout that it is enforced exactly. 
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where hk on the left hand side solves the mode equation Eq. ( |TTD exactly. Since for a 
generic finite rj* and fixed b, the recursion generated by Eq. eventually increases 

without bound in general, the recursion relation generates at best an asymptotic ex- 
pansion in the limit that the higher than zeroth adiabatic order terms go to zero. In 
particular, an infinite adiabatic order vacuum usually cannot be generated at a "nonsin- 
gular" 77*. 

Now, we examine how different boundary conditions (different adiabatic order vacua) 
give rise to different asymptotic behaviors as b — > 00. First let us restrict our attention 
to the case where a(t]) is C°° in the domain of interest. Since b 2 is the coefficient of a 2 
term inside w%, we see that a sufficient condition for the higher than zeroth adiabatic 
order terms to go to zero for large b is {d v a / drf) / a u+1 < 00 for all rj in the domain and 
any finite natural number v. Hence, we will assume this to be true and use the adiabatic 
expansion to determine the asymptotic power dependence of rix as 1/b — > 0. 

The key is that the recursion Eq. (A3) can be used as a generator of an asymptotic 



expansion of the exact solution in the limit that the higher than 0th adiabatic order terms 
tend to zero. One can easily show that this map has the property if A\Wk]/Wk ~ 
f + 0{l/b a ) + 0(l/b a+ ») with fi > 1, then A[A[w( n) }}/A[wt n) ] ~ f + 0{l/b a+2 ) + 
(9(l/6 a+2+M ) where ~ represents the asymptotic limit that b — > 00. Since wjf 1 = wj, + 
(9(1/5), we arrive at an useful property 

Wi n) = Wi n ~ 2) + 0{l/b n ~ l ), (A7) 

which shows how each successive approximation generates corrections of only increasingly 

higher order in 1/b. Thus, 

A/2-1 

W k {rj) ~ wM + J2 (Wf ^ - Wt ] ) + 0(l/b A+1 ) (A8) 

n=0 

is an asymptotic expansion of the solution to Eq. ( |A"2|) with the boundary condition 



WW) = wfV) + h( V *), = w' k (A \ V *) + h'( V *) (A9) 
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where h{rj) ~ 0{l/b A+1 ). Let us call this solution f k . Note that f k satisfies a boundary 
condition that differs from one implied by Eq. (|A6|) by 0(l/b A+1 ). 

To check this is the asymptotic expansion for the solution satisfying a different bound- 
ary condition (i.e. the one specified by Eq. (|A6p), one can now perturb about f k by writing 
Wfc (v) = fk(v) + u k{v) where the superscript on W k corresponds to the time at which 
the boundary condition Eq. ( |A6|) is imposed and by using this in Eq. (A.2) to obtain a 



differential equation linear in u k . One then finds that the sourced solution contributes 
only 0{l/b A+1 ) to u k (rj). Hence, if the initial data on u k {ri) is of the order of 0(l/b A+1 ), 
then the behavior of u k (rj) as b — > oo is 0(l/b A+1 ). In particular, if we have the boundary 
condition Eq. O instead of Eq. ([£]§, then we find u k (rj*) ~ 0{l/b A+l ) and 

WfW-W^^ + Oil/b^ 1 ) (A10) 

where the 0(l/b A+1 ) vanishes at r\ = r/*. This is of course what we would naively expect. 

We can now see how rix will depend asymptotically on b. Suppose the vacuum in the 
past is defined at rj = 770 with nth adiabatic order boundary condition and the vacuum 
today is defined at r] — r\\ with pth adiabatic order boundary condition. Carrying out 
the Bogoliubov transformation with the solution written in the form Eq. ( |Al|) , we find 



1 11/ w'J>° w h m ' 2 



(All) 



where the right hand side can be evaluated at any 77. In light of Eq. (|A10|) , if we substitute 
W? = Wi n) + 0{l/b n+1 ) and Wj] 1 = Wjf ] + C(1/6 P+1 ), then Wj! - = 0{l/b r+1 ) 
where r = Mm(n,p) and w' k m /W%° - W' k m jWf = 0(l/b r+2 ). Now, since 



oc 



n x oc / \l3 k \ 2 k z dk } (A12) 







after making a change of variable from k to y through k = yb, we obtain the result in 
Eq. (0). 
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If within the domain of interest there is one discontinuity of the first kind (left and 
right hand limits exist but are unequal) in (d q a/dr] q )/a q+1 at rj = r\ d for q = s where 
—2 < s — 2 < r = Min(p, n), and there are no discontinuities for q < s, then Eq. 
( |A11| ) will receive leading contributions at the discontinuity. Note that the asymptotic 
expansion is valid in each "continuous" region because the discontinuity is of the first 
kind. Hence, with similar considerations as with the smooth case above, we can obtain 
Eq. (|15|). However, unlike in the continuous case, the asymptotic expansion can be used 



to evaluate Eq. ( |A11|) only at r\ — r\ d because the asymptotic expansion cannot be 
extended beyond each of the continuous regions. If s is even, then W k ° (rj d ) — Wj] 1 (r] d ) ~ 
0{l/b s ~ l ) will give the leading contribution in Eq. ( |A11| ) because W"( s ~ 2 ^ /W^' 2 ^ is 
discontinuous. If s is odd, then the leading contribution to Eq. ( All|) will come from the 
difference W' k V0 {T] d )/W^{r] d )-W' k m {r] d )/W^{r] d ) ~ 0{l/b s ^) because (W"^ /W^)' 
is discontinuous. 
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